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A  molecular  dynamics  (MD)  simulation  is  used  to  reveal  the  grain  boundary  effect  on  the  ionic  transport 
of  yttria-stabilized  zirconia  (YSZ).  The  oxygen  ion  displacements  and  diffusivities  of  the  ideal  and  grain 
boundary-inserted  YSZ  models  are  analyzed  at  elevated  temperatures.  An  optimized  Y2O3  concentration 
within  YSZ  for  the  best  ionic  conductivity  is  achieved  by  balancing  the  trade-off  between  the  increased 
vacancies  and  the  decreased  accessible  free  space.  The  mass  transfer  resistance  of  the  grain  boundary 
in  YSZ  can  be  more  easily  found  at  higher  temperatures  by  observing  the  oxygen  ion  diffusivities  or 
traveling  trajectories.  At  lower  temperatures,  the  grain  interior  and  the  grain  boundary  control  the  ionic 
transport.  In  contrast,  the  grain  boundary  effect  on  the  diffusion  barrier  is  gradually  eliminated  at  elevated 
temperatures.  The  modeled  results  in  this  work  agree  well  with  previous  experimental  data. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Yttria-stabilized  zirconia  (YSZ)  solid  electrolytes  and  materials 
derived  from  them  have  widely  been  used  for  solid  oxide  fuel  cell 
applications  because  of  their  excellent  durability  for  long-term  use. 
In  order  to  improve  the  balance  between  performance  and  durabil¬ 
ity,  previous  studies  have  attempted  to  design  this  material  from 
the  perspective  of  the  material  composition  of  zirconia  [1-14].  The 
grain  boundary  (GB)  in  the  resultant  solid  electrolyte  exhibits  an 
apparently  higher  resistance  than  the  grain  interior  (GI)  [15],  which 
could  be  attributed  to  the  aggregation  of  complex  defects  near  the 
grain-grain  interface  regions  [16].  Therefore,  the  presence  of  GBs 
lowers  the  overall  ionic  conductivity  because  they  are  intrinsically 
unfavorable  to  ion  transport.  As  a  result,  previous  studies  have  ana¬ 
lyzed  ionic  transport  in  the  GIs  and  GBs  of  zirconia-based  solid 
electrolytes  [17-19];  however,  some  studies  have  suggested  that 
the  interfacial  area  between  two  grains  may  form  a  channel  for  ionic 
transport  with  the  proper  structural  design  [20,21  ],  which  also  indi¬ 
cates  the  importance  of  the  properties  of  GBs  for  ionic  transport.  As 
a  result,  GBs  also  play  a  vital  role  in  the  ionic  conduction  of  the  solid 
electrolyte  at  an  intermediate  temperature  for  SOFC  applications; 
however,  the  transport  mechanisms  of  oxygen  ions  across  the  GI 
and  GBs  are  not  easy  to  investigate  at  an  atomic  scale  using  exper¬ 
imental  methods.  Molecular  simulation  techniques  have  proven  to 
be  a  promising  tool  in  the  characterization  of  organic,  inorganic, 
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and  hybrid  materials  at  a  microscopic  level  [22-30].  The  molecular 
dynamics  (MD)  technique  can  be  used  to  successfully  describe  how 
doped  species  alter  the  ionic  conductivities  of  zirconia-based  solid 
electrolytes.  The  discussed  factors  include  the  size  of  the  cations, 
the  barrier  edge,  the  defect  complex,  interionic  interactions,  and 
the  associated  external  forces  [29-34].  Additionally,  phase  compo¬ 
sition  and  structural  stability  were  analyzed  using  MD  simulations, 
and  the  results  were  observed  to  agree  well  with  experimental 
results  [35,36].  Some  MD  studies  have  constructed  molecular  mod¬ 
els  in  order  to  analyze  the  ion  displacement  phenomenon  of  the 
GBs  in  YSZ  at  temperatures  greater  than  1273  K,  which  suggests 
that  GBs  typically  block  oxygen  ion  transport  [37,38];  however, 
the  effect  of  GBs  on  ion  transport  in  the  intermediate  temperature 
region  on  a  microscopic  scale  is  not  fully  understood.  Thus,  this 
work  constructed  a  molecular  model  in  order  to  study  the  effect  of 
GBs  on  ionic  transport  mechanisms  in  YSZ.  We  developed  a  molec¬ 
ular  model  of  a  YSZ  structure  with  several  GBs,  such  that  the  ionic 
diffusion  mechanism  therein  was  analyzed  and  compared  with  an 
ideal  YSZ  model  (without  GBs).  In  this  work,  all  the  simulations 
were  processed  using  the  Materials  Studio  package  from  Accelrys. 
The  details  of  the  model  construction  and  the  physical  property 
analyses  are  described  in  the  following  section. 

2.  Theoretical  method 

2.1.  Model  construction 

The  zirconia-based  molecular  model  was  constructed  using  the 
cubic  zirconia  unit,  as  shown  in  Fig.  1(a).  In  this  work,  a  molecular 
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Fig.  1.  Schematic  expression  of  the  (a)  grain-free,  (b)  grain  boundary-inserted  molecular  models,  and  (c)  oxygen  ion  moving  trajectory  types  in  the  sub-lattice  of  cubic 
zirconia  structure. 


model  consisting  of  64  unit  cells  was  built  to  be  the  basis  of  the  YSZ 
model.  Subsequently,  different  amounts  of  Y203  were  added  to  our 
basis  model  to  replace  the  original  Zr02  units  so  as  to  form  the  YSZ 
structures  with  different  Y203  concentrations.  Specific  quantities 
of  vacancies  were  produced  using  the  following  process: 


YOi.5^Y"zr  +  |Oo  +  iv£,  (1) 


The  positions  of  the  Zr  ions  that  were  replaced  by  the  Y  ions 
were  randomly  assigned  in  the  YSZ  structures  with  different  dop¬ 
ing  concentrations.  Vacancies  were  formed  by  randomly  removing 
specific  amounts  of  oxygen  ions  in  the  system.  Five  independent 
molecular  models  were  constructed  and  were  adopted  in  order  to 
calculate  the  average  physical  properties. 

Six  types  of  YSZ  models  with  different  Y203  concentrations  were 
constructed  in  order  to  optimize  the  doping  concentration.  Table  1 
summarizes  the  model  parameters  in  detail. 
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Table  1 

Parameters  of  the  YSZ  molecular  models. 


Model 

Y203  (mol%a) 

Zr+4/Zr02  (no.) 

Y+3/Y203  (no.) 

O-2  (no.) 

V3  (no.) 

06-YSZ 

6.22 

226/226 

30/15 

497 

15 

07-YSZ 

7.11 

222/222 

34/17 

495 

17 

08-YSZ 

8.02 

218/218 

38/19 

493 

19 

09-YSZ 

8.94 

214/214 

42/21 

491 

21 

10-YSZ 

10.34 

208/208 

48/24 

488 

24 

a  Mol%  ofY203  =  100  x  (mole  no.  of  Y203)/(mole  no.  ofY203  +  mole  no.  of  Zr02). 


Later,  in  order  to  simulate  the  effect  of  grain  boundaries  on  ionic 
transport  in  the  YSZ  model,  we  manually  shifted  part  of  the  YSZ 
model  structure  in  the  direction  of  the  x-axis  so  as  to  create  slipped 
dislocations  in  the  molecular  model.  The  slipping  distance  was  set 
to  0.4436  A,  which  effectively  reflected  the  grain  boundary  effect 
on  the  YSZ  model  in  this  work.  Four  subregions  were  divided  by 
rearranging  the  ion  positions  in  the  space,  as  shown  in  Fig.  1(b).  It 
was  assumed  that  these  space  rearrangements  in  a  cubic  structure 
model  expressed  the  characteristics  of  a  polycrystalline  YSZ  struc¬ 
ture.  In  addition,  in  order  to  maintain  the  ‘inserted-grain  boundary 
structure’  for  the  duration  of  the  MD  simulation,  several  cations 
were  fixed  in  the  YSZ  model  during  the  calculation.  All  the  YSZ 
models  were  processed  using  a  500-ps  MD  duration  from  873  K  to 
1473  K  under  an  NVT  ensemble  (a  fixed  atom  number,  cell  volume, 
and  system  temperature). 


Fig.  2.  Ionic  conductivity  values  from  this  simulation  and  experimental  work  [35] 
at  1273  K. 


2.3.1.  Mean-squared  displacement  (MSD )  and  self-diffusivity 

The  diffusion  mechanism  and  self-diffusivity  of  oxygen  ions  in 
the  YSZ  system  were  analyzed  using  mean-squared  displacement 
(MSD),  which  can  be  defined  using  the  Einstein  relationship: 


2.2.  Potential  function 


The  Born-Mayer-Buckingham  potential  was  adopted  in  order  to 
illustrate  the  metal  oxide  structure  and  ionic  transport  behaviors 
of  the  binary  or  mixed  oxide  doped  zirconia  system,  as  follows: 


J(r,j)  =  A j-exp 


c  +  ^Coulomb? 

r ?. 
y 


(2) 


where  Ay,  Cy ,  and  py  are  the  coefficients  that  are  used  in  the  poten¬ 
tial  equation.  In  Eq.  (2),  the  interaction  potential  between  two 
atoms,  i  and  j,  including  the  van  der  Waals  (first  two  terms)  and 
Coulombic  (third  term)  energies,  are  described.  In  addition,  in  Eq. 
(2),  r y  represents  the  distance  between  ions  i  and  j,  and  ^coulomb 
represents  the  Coulombic  term.  The  associated  parameters  for  the 
theoretical  calculations  in  this  work  were  validated  using  empirical 
fittings  from  previous  studies  that  have  investigated  zirconia-based 
structures  [39,40],  as  shown  in  Table  2. 


2.3.  Physical  property  analysis 


N 

MSD  (t)=L^  <[r,( to  +  t)  -  r,(t0)]2>  =  B  +  6D  x  t  (3) 

1=1 


where  N  is  the  total  number  of  atoms;  r,-(to  +  t)  and  r,-(t0)  are  the 
positions  at  time  t0  +  t  and  time  t0,  respectively;  B  is  a  constant; 
and  D  is  the  self-diffusivity.  The  self-diffusivity  can  be  obtained  as 
the  slope  of  the  MSD  curve  versus  the  elapsed  time.  The  diffusivity 
ratios  of  oxygen  ions  between  the  M-07-YSZ  and  P-07-YSZ  models 
can  be  estimated  using  the  follow  equation: 


Diffusivity  ratio  = 


Oxygen  ion  diffusivity  of  P-07-YSZ 
Oxygen  ion  diffusivityof  M-07-YSZ  ’ 


(4) 


2.3.2.  Ionic  conductivity 

The  ionic  conductivity  can  be  estimated  from  the  diffusion 
coefficient  of  oxygen  ions  using  the  Nernst-Einstein  relationship 
expressed  by  the  following  equation: 


Nq2D 


(5) 


In  our  simulation,  five  independent  molecular  models  were  con¬ 
structed  in  order  to  analyze  physical  properties  and  obtain  average 
calculated  values.  The  reproducibility  was  also  confirmed  in  order 
to  validate  that  this  was  a  reasonable  simulation  procedure. 


where  a  is  the  ionic  conductivity,  N  is  the  number  of  oxygen  ions 
in  the  calculation  model,  q  is  the  electrical  charge  of  the  ions,  D 
is  the  diffusion  coefficient,  /<b  is  the  Boltzmann  constant,  T  is  the 
operation  temperature,  and  HR  is  the  Elaven  coefficient,  which  was 
set  at  0.65  in  this  study  [41  ]. 


Table  2 

The  potential  parameters  that  were  used  in  the  Born-Mayer-Buckingham  potential 
in  this  study. 


Ion  pair 

Ay  (kcalmol-1) 

Pi/ (A) 

Cy  (kcalA6  mop1) 

Reference 

Zr-Zr 

0 

1 

0 

Zr-Y 

0 

1 

0 

Y-Y 

0 

1 

0 

Zr-0 

22717.937 

0.376 

0 

[40] 

Y-0 

30995.940 

0.349 

0 

[39] 

0-0 

524480.535 

0.149 

642.704 

[40] 

2.3.3.  Ion  movement  distribution 

During  the  MD  simulation,  the  actual  ionic  moving  distance  of 
each  oxygen  ion  can  be  estimated  using  its  initial  and  final  coordi¬ 
nates.  The  statistics  of  each  moving  distance  and  the  corresponding 
quantity  were  collected  to  express  the  distribution  of  ion  move¬ 
ment.  The  moving  distances  of  ions  can  be  compared  with  the 
oxygen-oxygen  distances  in  the  cubic  zirconia  unit,  which  can 
then  be  categorized  as  one  of  several  different  movement  types, 
as  shown  in  Fig.  1(c). 
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Fig.  3.  Accessible  free  space  and  vacancy  numbers  of  YSZ  models  with  various  Y2O3 
concentrations  at  1273  K. 


3.  Results  and  discussion 

3. 1.  Oxygen  ion  transport  in  ideal  YSZ  model 

Fig.  2  shows  the  oxygen  ionic  conductivities  of  ideal  YSZ  models 
with  different  doped  Y203  concentrations.  In  Fig.  2,  it  can  be  clearly 
seen  that  the  07-YSZ  revealed  the  highest  ionic  conductivity  among 
all  the  YSZ  models  with  different  compositions.  In  the  YSZ  struc¬ 
ture,  increasing  the  doped  Y203  would  release  more  vacancies.  The 
increased  vacancies  from  adding  Y203  to  the  YSZ  structure  tended 
to  form  the  pathways  for  oxygen  ion  transport,  which  is  benefi¬ 
cial  for  improving  the  ionic  conductivity.  Flowever,  while  the  Y203 
concentration  is  higher  than  7  mol%,  there  is  a  paradox  between  the 
elevated  vacancy  number  and  the  ionic  conductivity,  as  shown  in 
Fig.  2.  The  optimized  Y203  concentration  of  YSZ  for  ionic  transport 
was  approximately  7  mol%  in  this  work.  From  the  past  studies,  it 
was  determined  that  the  excess  Y+3  in  the  YSZ  system  tended  to 
trap  the  vacancies,  which  lowered  the  vacancy  mobility.  Thus,  oxy¬ 
gen  ion  mobility  was  also  inhibited,  which  led  to  a  less  effective 
thermal  motion  [30,31  ].  In  addition,  some  studies  also  pointed  out 
that  the  Y-Y  edge  formed  in  the  YSZ  system  also  hampered  the  oxy¬ 
gen  ion  movement  and  then  decreased  the  ion  mobility  [29-31,42]. 
As  a  result,  an  optimized  Y203  concentration  for  the  YSZ  struc¬ 
ture  would  be  found.  Meanwhile,  we  also  found  that  the  07-YSZ 
model  provided  the  best  ionic  conductivity,  which  was  similar  to 
previous  experimental  results  [43],  as  shown  in  Fig.  2.  Thus,  the 
07-YSZ  model  was  adapted  for  the  analysis  of  the  impact  of  grain 
boundaries  on  ionic  transport  behavior  in  this  study.  Flowever,  we 
can  still  found  some  differences  between  the  simulated  and  exper¬ 
imental  ionic  conductivities  of  YSZ  at  various  Y203  concentrations. 
The  reasons  of  causing  this  error  in  our  simulation  work  might  be 
explained  by  the  following  descriptions.  Firstly,  the  dimension  of 
the  simulated  model  and  the  MD  duration  adopted  in  this  work 
are  less  than  that  of  the  real  system  due  to  the  limited  computa¬ 
tional  resources  presently.  Meanwhile,  some  intrinsic  errors  in  the 
calculation  process  might  still  be  produced  by  the  selected  poten¬ 
tial  function  and  parameters  optimized  from  the  semi-empirical 
equations.  Flere,  the  comparable  quantitative  results  between  the 
simulation  and  experimental  works  can  be  observed,  which  sug¬ 
gests  that  the  MD  simulation  can  be  a  useful  for  the  prediction  of 
ionic  transport  ability  of  zirconia-based  solid  electrolytes. 

Meanwhile,  we  also  discussed  the  optimized  Y203  concentra¬ 
tion  in  YSZ  for  the  best  ionic  conductivity  from  the  perspective 
of  free  volume.  Fig.  3  illustrates  the  number  of  vacancies  and  the 
accessible  free  space  (AFS)  corresponding  to  the  Y203  concentra¬ 


Fig.  4.  MSD  diagrams  of  zirconium,  yttrium,  and  oxygen  ions  in  07-YSZ  at  1273  K. 

tions  of  YSZ  models.  The  AFS  value  can  be  obtained  by  probing  the 
available  space  for  a  hard,  spherical  particle  passing  through  a  YSZ 
model,  where  the  radius  of  the  probe  was  set  as  that  of  the  oxygen 
ion.  A  higher  AFS  value  indicates  more  and  larger  internal  spaces 
or  a  higher  connectivity  among  the  vacancies  in  the  YSZ  structure, 
which  might  improve  ionic  transport.  Flowever,  it  should  also  be 
noted  that  the  AFS  decreased  with  the  increase  of  Y203  in  the  YSZ 
models.  This  reduction  of  AFS  with  the  increase  of  vacancies  sug¬ 
gested  that  the  connectivity  of  vacancies  might  be  reduced,  or  the 
pathway  of  oxygen  ions  might  be  hindered  by  raising  the  Y203  con¬ 
centration.  Moreover,  the  distribution  of  vacancies  in  higher  YSZ 
concentration  systems  might  not  be  suitable  for  the  formation  of 
continuous  ion  motion,  which  is  not  favorable  when  trying  to  pro¬ 
duce  better  ionic  transport.  In  contrast,  the  increased  Y203  content 
would  also  enlarge  the  lattice  constant  of  the  YSZ  structure  because 
of  the  larger  radius  of  the  yttrium  ion,  which  led  to  the  expansion 
of  the  molecular  model.  The  expansion  of  the  YSZ  structure  might 
also  increase  the  free  space  in  the  system.  Thus,  we  checked  the 
lattice  parameters  of  6  mol%  and  10  mol%  YSZ  from  the  past  theo¬ 
retical  and  experimental  works  [31  ].  The  ratio  of  lattice  parameter 
of  10  and  6mol%  YSZ  is  approximately  1.0019.  This  extent  of  free 
space  expansion  would  not  alter  the  results  of  the  AFS  analysis  in 
this  work,  as  shown  in  Fig.  3. 

In  addition,  the  displacements  of  zirconium  and  yttrium  ions  in 
the  07-YSZ  model  were  much  smaller  than  those  of  oxygen  ions,  as 
illustrated  in  Fig.  4.  During  the  500-ps  MD  simulation,  the  overall 
displacements  of  zirconium  and  yttrium  ions  were  both  approxi¬ 
mately  0.06  A,  whereas  those  of  oxygen  ions  were  much  higher  at 
1 0  A.  The  cations  may  have  been  nearly  fixed  or  may  have  fluctuated 
around  the  lattice  sites.  Thus,  it  was  assumed  that  partially  fixed 
cations  did  not  affect  the  transport  of  oxygen  ions  in  the  YSZ  model. 
In  this  work,  the  grain  boundary  model  is  referred  to  as  ‘P-07-YSZ’, 
whereas  the  ideal  YSZ  model  is  referred  to  as  ‘M-07-YSZ’. 

3.2.  GB  effect  on  oxygen  ion  MSD  and  movement 

To  compare  the  ion  mobility  in  YSZ  structures  between  the  ideal 
and  GB-inserted  models,  their  MSD  diagrams  were  analyzed,  as 
illustrated  in  Fig.  5.  In  Fig.  5,  it  can  be  clearly  seen  that  the  P-07- 
YSZ  model  showed  lower  MSD  slopes  and  displacements  during  the 
500-ps  calculation.  Because  the  P-07-YSZ  model  was  also  composed 
of  pure  cubic-phase  units,  the  structure  of  the  P-07-YSZ  model  was 
no  different  from  that  of  M-07-YSZ,  except  for  the  presence  of  a 
“grain  boundary.”  Therefore,  this  lower  ion  mobility  of  P-07-YSZ 
might  be  caused  by  the  migration  barrier  from  the  inserted  GBs. 
The  slipped  dislocations  inside  P-07-YSZ  formed  the  ion  transport 
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Fig.  5.  MSD  diagrams  of  the  (a)  M-07-YSZ  and  (b)  P-07-YSZ  models  at  different 
temperatures. 


Oxygen  Ions 


b 


Oxygen  Ions 

Fig.  6.  Oxygen  ion  displacements  during  500-ps  MD  simulations  using  the  (a)  M- 
07-YSZ  and  (b)  P-07-YSZ  models  at  1273  K. 

GBs.  Even  though  vacancies  or  mobile  oxygen  ions  were  formed 
near  the  GBs  during  the  model  creation  process,  the  higher  mass 
resistances  of  the  GBs  impeded  their  ability  to  travel  across  the 
boundary.  Therefore,  the  percentage  of  immobile  oxygen  ions  was 
higher  because  of  the  GBs.  In  contrast,  the  M-07-YSZ  model  exhib- 


resistance,  which  hindered  partial  ionic  hopping  from  one  site  to 
another.  Thus,  the  concept  of  inserted  grain  boundary  can  be  vali¬ 
dated  here  via  the  MSD  diagram  analysis.  In  order  to  compare  the 
transport  behaviors  of  oxygen  ions  in  the  P-07-YSZ  model  with 
those  in  the  M-07-YSZ  model,  we  estimated  the  absolute  values  of 
displacements  of  all  oxygen  ions  in  both  models.  Fig.  6  illustrates 
the  oxygen  ion  displacement  values  of  the  M-07-YSZ  and  P-07-YSZ 
models.  These  displacement  values  can  be  roughly  classified  into 
three  categories:  below  1  A,  between  2  and  6  A,  and  beyond  6  A. 
In  the  region  between  2  and  6  A,  it  can  be  observed  that  M-07-YSZ 
had  a  higher  value  than  that  of  P-07-YSZ.  M-07-YSZ  contained  more 
displacements  with  high  values  in  comparison  with  those  of  P-07- 
YSZ,  which  indicates  its  higher  transport  ability  for  oxygen  ions. 
Moreover,  a  number  of  oxygen  ions  with  long-distance  movement 
(larger  than  8  A)  can  be  observed  in  M-07-YSZ  (Fig.  6(a)),  whereas 
no  such  oxygen  ions  can  be  observed  in  P-07-YSZ  (Fig.  6(b)).  Thus, 
it  can  be  inferred  that  the  GBs  resisted  ionic  transport  and  inhibited 
continuous  ionic  movements. 

Fig.  7  depicts  the  displacement  distribution  of  oxygen  ions  in 
the  M-07-YSZ  and  P-07-YSZ  models.  In  Fig.  7,  the  P-07-YSZ  model 
exhibits  a  higher  percentage  of  displacement  values  that  are  located 
in  the  region  of  0-1 A  in  comparison  with  the  M-07-YSZ  model. 
More  specifically,  a  higher  number  of  ion  vibrations  or  ion  fluctu¬ 
ations  were  observed  in  the  YSZ  model  with  GBs.  This  higher  ion 
immobility  can  be  attributed  to  the  barrier  that  is  formed  by  the 


Fig.  7.  Comparison  of  the  oxygen  ion  displacement  profiles  of  the  M-07-YSZ  and 
P-07-YSZ  models  during  500-ps  MD  simulations  at  1273 1<. 
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Temperature  (K) 

Fig.  8.  Oxygen  ion  diffusivities  of  the  (a)  M-07-YSZ  and  (b)  P-07-YSZ  models  at 
different  temperatures. 

ited  more  high-displacement  ions  (a  displacement  value  of  greater 
than  3  A)  than  the  P-07-YSZ  model,  which  indicates  better  ionic 
transport  in  the  ideal  model.  In  comparison  with  the  ideal  cubic 
YSZ  model,  the  model  with  the  GBs  suggested  a  method  of  oxy¬ 
gen  ion  transport  exists  in  the  YSZ  structure,  with  results  that  were 
similar  to  previously  obtained  experimental  results. 

3.3.  GB  effect  on  oxygen  ion  diffusivity 

The  diffusivities  of  oxygen  ions  in  the  M-07-YSZ  and  P-07-YSZ 
models  were  also  estimated  in  this  study,  as  illustrated  in  Fig.  8. 
In  Fig.  8,  the  oxygen  ion  diffusivities  of  the  two  models  were 
improved  with  increases  in  temperature.  Moreover,  higher  tem¬ 
peratures  were  observed  to  increase  the  kinetic  energy,  which  also 
increases  ion  velocity  and  directs  ion  movement.  Furthermore,  the 
differences  in  the  oxygen  ion  diffusivities  between  the  M-07-YSZ 
and  P-07-YSZ  models  were  observed  to  increase  with  increased 
temperature,  which  indicates  that  the  effect  of  the  GBs  on  oxygen 
ion  transport  can  be  more  easily  observed  at  elevated  temper¬ 
atures.  The  enhanced  differences  in  the  oxygen  ion  diffusivities 
that  accompanied  elevated  temperatures  may  be  explained  by  the 
following  phenomenon.  In  the  case  of  a  fluorite-structured  solid 
electrolyte,  ionic  conduction  is  induced  when  the  material  is  heated 
to  1073  K  (Ce02-based)  or  1273  K  (Zr02-based).  In  other  words,  at 


Fig.  9.  Oxygen  ion  diffusivity  ratios  of  the  P-07-YSZ  to  M-07-YSZ  models  at  different 
temperatures. 


lower  temperatures,  ion  migrations  might  not  be  clearly  observ¬ 
able.  In  this  situation,  both  the  M-07-YSZ  and  P-07-YSZ  models 
exhibited  lower  ion  diffusivities.  The  effect  of  the  GBs  on  diffusivity 
is  therefore  not  apparent  in  the  P-07-YSZ  model.  At  higher  tem¬ 
peratures,  the  oxygen  ions  in  the  P-07-YSZ  model  demonstrated 
enhanced  long-distance  transport.  The  probability  of  ion  interac¬ 
tion  with  the  GBs  also  increased.  Therefore,  intergrain  resistance 
was  dominant  in  this  scenario. 

In  contrast,  in  the  M-07-YSZ  model,  oxygen  ion  transport  was 
primarily  impacted  by  GI  resistance,  and  a  higher  energy  (due  to  a 
higher  temperature)  would  permit  the  ions  to  overcome  the  intra¬ 
grain  barrier,  regardless  of  temperature.  Therefore,  the  differences 
between  the  ideal  and  the  GBs-inserted  model  were  substantial 
with  increases  in  temperature.  In  addition,  we  observed  a  break 
point  at  approximately  1 273 1<  by  using  a  regression  of  the  diffusiv¬ 
ity  values  (M-07-YSZ)  in  Fig.  8,  which  approached  the  temperature 
that  is  required  to  induce  a  high  diffusivity.  A  similar  result  was 
observed  in  the  P-07-YSZ  model.  The  obvious  increases  in  diffusiv¬ 
ity  indicate  that  the  impact  of  GBs  on  ion  migration  might  be  easier 
to  observe  at  higher  temperatures. 

Unfortunately,  it  is  difficult  to  discuss  the  extent  of  the  influ¬ 
ence  of  the  GBs  on  ion  migration  by  directly  observing  differences 
in  diffusivity.  A  direct  comparison  of  diffusivity  differences  may  not 
consider  how  the  ionic  transport  capability  of  the  P-07-YSZ  model 
differs  from  that  of  the  M-07-YSZ  model  based  on  the  performance 
of  the  M-07-YSZ  model.  Therefore,  in  order  to  explain  the  asso¬ 
ciated  ion  transport  mechanism,  we  calculated  the  ratios  of  the 


Temperature  (K) 


Temperature  (K) 


Fig.  10.  Differences  in  the  oxygen  ion  diffusivity  of  the  (a)  M-07-YSZ  and  (b)  P-07- 
YSZ  models  at  each  temperature  interval  from  873  K  to  1473  K. 
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Suppositional  Boundary 


M-07-YSZ,  t  =  0  ps 


Grain  Boundary 


P-07-YSZ,  t  =  0  ps 


M-07-YSZ,T  =  773  K,  t  =  500  ps 


M-07-YSZ,T  =  873  K,  t  =  500  ps 


M-07-YSZ,T  =  973  K,  t  =  500 ps 


P-07-YSZ,  T  =  773  K,  t  =  500  ps 


P-07-YSZ,  T  =  873  K.  t  =  500  ps 


P-07-YSZ,T  =  973  K,t  =  500  ps 


M-07-Y5Z,  T  =  1073  K,  t  =  500  p$ 


M-07-YSZ,  T  =  1173  K,  t  =  500  ps 


Fig.  11.  Traveling  trajectories  of  oxygen  ions  in  M-07-YSZ  and  P-07-YSZ  models  from  773 1<  to  1473 1<. 


oxygen  ion  diffusivities  of  the  two  investigated  models  at  differ¬ 
ent  temperatures.  Fig.  9  illustrates  the  diffusivity  ratios  of  oxygen 
ions  in  the  M-07-YSZ  and  P-07-YSZ  models.  The  diffusivity  ratio 
indicates  the  variation  in  ion  diffusivity  between  the  ideal  and  GB- 
inserted  models.  When  this  ratio  approaches  unity,  it  implies  that 


the  GBs  have  no  effect  on  ionic  transport  (no  difference  between 
the  M-07-YSZ  and  P-07-YSZ  models);  however,  a  decreased  ratio 
suggests  that  the  GBs  in  the  P-07-YSZ  model  do  affect  ionic  trans¬ 
port.  In  Fig.  9,  the  diffusivity  ratio  was  first  decreased  (in  Region 
I)  and  then  increased  (in  Region  II)  by  elevating  the  operational 
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P-07-YSZ,T  =  1273  K,  t  =  500  ps 


P-Q7-YSZ,T  =  1373  K,  t  =  500  ps 


Fig.  11.  ( continued ). 


temperature.  This  change  in  diffusivity  ratio  might  be  explained  by 
the  following  phenomenon.  In  real  YSZ-based  materials,  the  energy 
barrier  to  ion  transport  can  be  classified  into  two  categories:  disso¬ 
ciation  energy  and  migration  energy  [36].  The  dissociation  energy 
indicates  that  the  ions/vacancies  that  are  bound  to  the  surrounding 
ions  trap  the  ions/vacancies.  The  ion  can  leave  the  trapped  center 
when  it  acquires  sufficient  energy  such  that  it  can  break  away  from 
its  bound  state.  The  migration  energy  indicates  that  the  energy  bar¬ 
rier  that  is  formed  by  the  resistance  of  ions  that  travel  between 
the  two  sites  might  be  formed  by  an  intrinsic  transport  resistance, 
geometric  structural  factors  (such  as  the  presence  of  a  GB),  unfa¬ 
vorable  phase  structures,  or  some  defects  in  the  system.  In  Region 
I  (lower  temperature),  the  defect  complexes  consist  of  vacancies 
and  doped  cations  with  high  binding  energies  that  lower  vacancy 
mobility.  These  defect  complexes  tend  to  accumulate  near  the  GBs, 
which  also  increase  the  resistance.  Meanwhile,  the  energy  barrier  of 
oxygen  ions  moving  from  one  site  to  another  would  also  be  higher. 
Increasing  the  operational  temperature  increases  the  number  of 
ions  that  would  interact  with  GBs,  which  reflects  the  influence  of 
GBs  on  diffusivity.  Thus,  the  influence  of  GBs  on  ion  transport  in  the 
P-07-YSZ  model  gradually  increased  with  increasing  temperature 
in  Region  I,  which  was  reflected  by  the  reduced  diffusivity  ratios.  In 
Region  II,  the  elevated  temperature  provides  more  kinetic  energy 
to  the  vacancies,  which  may  decrease  the  frequency  of  defect  com¬ 
plex  formation.  Meanwhile,  ions  with  higher  velocities  were  able 


to  overcome  the  GB-induced  geometric  resistance.  Therefore,  as 
the  system  was  heated  to  specific  temperatures,  the  associated 
mass  transport  resistances  may  have  changed.  The  GB-induced 
oxygen  ion  migrational  energy  barrier  could  be  gradually  reduced 
as  a  result  of  (a)  the  ‘decomposition’  of  defect  complexes  in  the 
system  (including  those  near  the  GBs)  and  (b)  increased  ion  mobil¬ 
ity  at  higher  temperatures.  In  other  words,  the  impact  of  GBs  on 
ion  transport  gradually  decreased  in  the  P-07-YSZ  model  at  ele¬ 
vated  temperatures  in  Region  II.  As  a  result,  the  diffusivity  ratio 
tended  to  increase  with  increasing  temperature  in  Region  II.  This 
phenomenon  of  a  reduced  impact  of  GB  on  ion  transport  at  high 
temperatures  also  agrees  well  with  previous  experimental  results 
[4,10,16]. 

In  addition,  we  also  estimated  the  difference  in  ion 
diffusivities  between  the  two  specific  temperatures,  e.g., 
ADs73-973  =  Dg73  —  D873,  wherein  Dg 73  and  are  the  diffu¬ 
sivities  of  oxygen  ions  at  873  K  and  973  K,  respectively,  as  shown 
in  Fig.  10.  As  shown  in  Fig.  10(b),  the  AD  values  tended  to  be 
maintained  or  slightly  reduced  at  low  temperatures  and  increased 
at  high  temperatures  in  the  P-07-YSZ  model.  This  tendency 
illustrates  that  the  improvement  in  ion  diffusivity  via  heating  at 
low  temperatures  was  limited  because  of  the  apparent  effect  of 
GBs.  Meanwhile,  the  AD  values  in  the  P-07-YSZ  model  increased 
at  higher  temperatures  because  the  mass  transport  barrier  was 
reduced  in  the  absence  of  GB  resistance;  however,  the  AD  values 
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in  the  M-07-YSZ  model  continuously  increased  throughout  the 
entire  investigated  temperature  range,  which  indicates  that  the 
ideal  YSZ  model  only  revealed  the  mass  transport  resistance  of 
the  GI.  Thus,  the  AD  values  were  monotonically  increased  by  the 
heating  treatment  and  did  not  indicate  the  impact  of  the  added 
GBs. 

3.4.  Oxygen  ion  traveling  trajectory 

We  also  observed  the  ion  traveling  trajectories  of  M-07-YSZ  and 
P-07-YSZ  models  during  the  500-ps  MD  duration  at  different  tem¬ 
peratures,  as  Fig.  11  shows.  In  Fig.  11,  we  selected  several  oxygen 
ions  in  both  YSZ  models.  The  initial  positions  of  these  ions  were 
set  as  the  same  coordinates  for  comparing  their  trajectories.  The 
first  row  shows  the  ion  positions  of  YSZ  models  at  the  initial  point 
of  the  MD  calculation,  i.e.,  t=0ps.  The  following  rows  show  the 
ion  positions  of  YSZ  models  after  500-ps  MD  calculations.  In  other 
words,  the  images  show  the  ion  positions  at  t=500ps  at  various 
temperatures.  In  Fig.  11,  it  was  found  that  the  ions  neither  left  the 
suppositional  boundary  (M-07-YSZ)  nor  the  grain  boundary  (P-07- 
YSZ)  while  the  temperature  was  below  973  K.  This  slight  oxygen 
movement  pointed  out  the  minute  ion  mobility  at  low  temperature. 
Furthermore,  the  ions  started  to  show  a  clear  migration  behavior 
from  1073  K  to  1274K.  The  more  drastic  ion  transport  was  elicited 
when  the  temperature  was  higher  than  1373  K.  Furthermore,  it 
should  also  be  noted  that  the  difference  in  ion  movements  between 
M-07-YSZ  and  P-07-YSZ  became  much  more  obvious  at  higher  tem¬ 
peratures,  especially  in  the  case  of  1273  K,  1373  K,  and  1473  K. 
The  higher  temperature  made  the  GB  effect  easier  to  observe  by 
tracing  the  ion  positions.  This  temperature-dependent  ionic  trans¬ 
port  agreed  with  the  results  of  the  diffusivity  analysis  that  were 
described  above  (see  Fig.  8). 

4.  Conclusions 

In  this  work,  a  molecular  YSZ  model  with  added  GBs  was  suc¬ 
cessfully  constructed  via  an  MD  simulation  technique.  An  increase 
in  the  Y203  concentration  in  YSZ  released  more  vacancies  but 
also  lowered  the  accessible  free  space  and  vacancy  mobility.  These 
opposite  tendencies  caused  by  increasing  the  Y203  content  led  to 
an  optimized  YSZ  concentration  at  7  mol%.  The  presence  of  GBs  in 
the  YSZ  structure  resulted  in  degraded  ion  transport  characteristics, 
which  increased  low-mobility  ions  and  decreased  continuous  ion 
motions  near  the  boundaries,  resulting  in  lower  ion  diffusivities. 
The  GB  effect  can  be  more  easily  observed  at  higher  temperatures. 
From  the  ion  diffusivity  analyses,  the  GI  and  GBs  dominate  ionic 
transport  at  lower  temperatures,  whereas  the  influence  of  GBs  was 
eliminated  at  higher  temperatures  because  the  migration  energy 
barrier  was  overcome.  The  GB-inserted  YSZ  model  can  effectively 
simulate  ionic  transport  through  the  GI  and  GBs,  thus  providing  a 
more  accurate  methodology  for  predicting  ionic  transport  behav¬ 
iors  in  comparison  with  that  achievable  using  a  conventional  ideal 
YSZ  model. 
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